\documentclass{sig-alternate}

\begin{document}

\title{L'Ecuyer y tiempo de vuelo del Enterprise}

\numberofauthors{2}

\author{
    \alignauthor
    Mui\~{n}a, Pablo\\
    \email{pmuina@alu.itba.edu.ar} \\
    \ \\
    Valdivieso, Ignacio\\
    \email{ivaldivi@alu.itba.edu.ar} \\
}

\maketitle

\begin{abstract}
En este art\'{i}culo se estudia el generador de n\'{u}meros pseudoaleatorios
de L'Ecuyer y se lo pone a prueba utilizando los tests de $\chi^{2}$ y de 
Kolmogorov-Smirnov. Luego se lo utiliza para generar n\'{u}meros pseudoaleatorios
con los que se estima el tiempo de vuelo medio del USS Enterprise mediante
simulaciones con el m\'{e}todo de Montecarlo.
\end{abstract}

\keywords{L'Ecuyer, $\chi^{2}$, Kolmogorov-Smirnov,
simulaci\'{o}n de Montecarlo, USS Enterprise}

\section{Introducci\'{o}n}\label{introduccion}

Los n\'{u}meros aleatorios surgen \'{u}nicamente de procesos naturales, los 
cuales no pueden ser recreados en un dispositivo determinista como una 
computadora. Sin embargo es posible generar secuencias de n\'{u}meros suficientemente
largas tal que no sea posible distinguir que los mismos no fueron generados
de manera aleatoria. Estos se conocen como n\'{u}meros pseudoaleatorios.

En la secci\'{o}n \ref{generador} se presenta el generador de L'Ecuyer, que 
utiliza dos LCG para generar las secuencias. Se analiza las secuencias que
genera el mismo, y se lo somete a tests estad\'{i}sticos para analizar la
distribuci\'{o}n de los n\'{u}meros generados. En la secci\'{o}n \ref{triangular}
se transforma la distribuci\'{o}n de la secuencia de n\'{u}meros obtenidos
por el generador de uniforme a triangular mediante la t\'{e}cnica de la
transformada inversa. En la secci\'{o}n \ref{simulacionpropulsor} se realiza
una estimaci\'{o}n del tiempo medio de vuelo de la nave USS Enterprise
utilizando el m\'{e}todo de Montecarlo.

\section{Generador de L'Ecuyer}\label{generador}

El generador de L'Ecuyer combina dos LCG para la generaci\'{o}n de n\'{u}meros 
pseudoaleatorios, de acuerdo al siguiente algoritmo:

$1$. Seleccionar una semilla $X_{1,0}$ en el rango $[1, 2147483562]$ para el 
LCG$1$ y $X_{2,0}$ en el rango $[1, 2147483398]$ para el LCG$2$.

$2$. Evaluar cada generador individual

\begin{equation}
\label{LCG1}
X_{1,n+1} = 40014 \ X_{1,n} \ mod \ 2147483563
\end{equation}

\begin{equation}
\label{LCG2}
X_{2,n+1} = 40692 \ X_{2,n} \ mod \ 2147483399
\end{equation}

$3$. Computar

\begin{equation}
\label{xn}
X_{n+1} = (X_{1,n+1} - X_{2,n+1}) \ mod \ 2147483562
\end{equation}

$4$. Computar

\begin{equation}
\label{un}
U_{n+1} = 
    \begin{cases}
    \frac{X_{n+1}}{2147483563}, X_{n+1} > 0\\
     \ \\
    \frac{2147483562}{2147483563}, X_{n+1} = 0\\
    \end{cases}
\end{equation}

$5$. Hacer $n = n + 1$ e ir a $2$.

A partir de este algoritmo se generan 10000 n\'umeros pseudoaleatorios, usando como semillas $X_{1,0}=3543$ y $X_{2,0}=2349$. Estos n\'umeros 
obtenidos se dividen en 10 intervalos de clase del mismo ancho. Las frecuencias obtenidas se pueden ver en el histograma de la figura \ref{fig:histogramalecuyer}, que muestra una salida uniforme. La decisi\'{o}n de elegir 10 intervalos de clase se basa en la proposici\'{o}n
de Nu\~{n}es (1985): la cantidad necesaria para que el histograma se vea lindo.

\begin{figure}[t]
\label{fig:histogramalecuyer}
\centering
\includegraphics[width=3.2in]{histoLecuyer}
\caption{Frecuencia de 10000 n\'umeros generados con el generador de L'Ecuyer divididos en 10 intervalos de clase}
\end{figure}


Se analizan gr\'{a}ficamente los numeros obtenidas con el algoritmo de L'Ecuyer para determinar la dependencia entre realizaciones. En la figura 2 se grafican las duplas $(U_i, U_{i+1})$, y en la figura 3 las ternas $(U_i, U_{i+1}, U_{i+2})$. Como se puede ver en ambas figuras, no se detecta
a simple vista que las duplas o las ternas se dispongan en forma de hiperplanos (rectas en el caso de las duplas y planos en el caso de las ternas), sino que se disponen de forma que hace que la salida parezca aleatoria.


\begin{figure}[t]
\label{fig:2d}
\centering
\includegraphics[width=3.2in]{duplas}
\caption{Duplas $(U_i, U_{i+1})$ de numeros generados por el algoritmo de L'Ecuyer}
\end{figure}

\begin{figure}[t]
\label{fig:3d}
\centering
\includegraphics[width=3.2in]{terna}
\caption{Ternas $(U_i, U_{i+1}, U_{i+2})$ de numeros generados por el algoritmo de L'Ecuyer}
\end{figure}

%\section{Probando el generador de L'Ecuyer}\label{pruebagenerador}

Para verificar si la secuencia generada se encuentra
distribuida de forma uniforme, se utiliza el test $\chi^{2}$ y el test de 
Kolmogorov-Smirnov.

\subsection{Test $\chi^{2}$}\label{testchicuadrado}

Se desea aceptar la hip\'{o}tesis 
$H_{0}:$ la secuencia de n\'{u}meros pseudoaleatorios generada utilizando el 
algoritmo de L'Ecuyer se encuentra uniformemente distribuida en el
intervalo $(0,1)$ con un nivel de confianza del $95\%$.

Por tratarse de una distribuci\'{o}n uniforme y al tomar 10 intervalos de clase,
la media esperada en cada uno es $E_{i} = 1000$. Entonces el estad\'{i}stico dado por la 
ecuaci\'{o}n:

\begin{equation}\label{estadistico}
\chi^{2}_{0} = \sum_{i=1}^{100} \frac{(O_{i} - E_{i})^{2}}{E_{i}}
\end{equation}

resulta aproximadamente $7.1980$. Como se toman 10 intervalos, la distribuci\'{o}n $\chi^{2}$
posee $9$ grados de libertad. Para $\alpha = 0.05$, se obtiene el valor cr\'{i}tico 
$\chi^{2}_{99,0.05} = 16.919$. Como 
$\chi^{2}_{99,0.05} > \chi^{2}_{0}$ se acepta $H_{0}$.

\subsection{Test Kolmogorov-Smirnov}\label{testks}

En el test de Kolmogorov-Smirnov se compara la distribuci\'{o}n emp\'{i}rica
obtenida a partir del generador con la distribuci\'{o}n te\'{o}rica de la variable aleatoria analizada.

Se utiliza la hip\'{o}tesis de que si se cumple que $D < D_{\alpha}$ los n\'{u}meros pseudoaleatorios generados
utilizando el algoritmo de L'Ecuyer se encuentran uniformemente distribuidos, donde $D$ es el
resultado de la prueba y $D_{\alpha}$ es el valor cr\'itico correspondiente a los par\'ametros de la prueba.

Realizando los calculos se obtiene $D = 0.1$, y para un nivel de significaci\'{o}n de $0.05$ se tiene 
$D_{\alpha} = 0.410$, por lo que se acepta la hip\'{o}tesis.

\section{Densidad triangular}\label{triangular}

En esta secci\'{o}n se utiliza la t\'{e}cnica de la transformada inversa, 
aplicada a la siguiente funci\'{o}n de densidad triangular:

\begin{equation}
\label{triangularEc}
f_{X}(x) = 
    \begin{cases}
    \frac{2(x-a)}{(b-a)(c-a)}, &  a \le x \le b\\
     \ \\
    \frac{2(c-x)}{(c-b)(c-a)}, & b < x \le c\\
     \ \\
     0, & \text{en otro lado}\\
    \end{cases}
\end{equation}

Esta t\'{e}cnica establece que dado un n\'{u}mero $u$ que es una realizaci\'{o}n
de $\mathcal{U} \sim \mathcal{U}[0,1]$, entonces $F^{-1}(u)$ es una
realizaci\'{o}n de una variable aleatoria $X$ con funci\'{o}n de distribuci\'{o}n
$F(x)$.

Integrando \ref{triangularEc} se obtiene la funci\'{o}n de 
distribuci\'{o}n $F_{X}(x)$:

\begin{equation}
\label{triangularIntegrada}
F_{X}(x) = 
    \begin{cases}
    0, & x < a\\
     \ \\
    \frac{x^{2} - 2ax + a^{2}}{(b-a)(c-a)}, &  a \le x \le b\\
     \ \\
    \frac{b-a}{c-a} + \frac{2cx - x^{2} - b(2c-b)}{(c-b)(c-a)}, & b < x \le c\\
     \ \\
     1, & x > c\\
    \end{cases}
\end{equation}

Calculando la inversa se llega a:

\begin{equation}
\label{vaconftriangular}
X = 
    \begin{cases}
    a + \sqrt{U(c-a)(b-a)}, 0 \leq U \leq \frac{b-a}{c-a} \\
    \ \\
    c - \sqrt{(c-b)^{2} - U(c-b)(c-a) + (b-a)(c-b)}, \\ \frac{b-a}{c-a} < U \leq 1 \\
    \end{cases}
\end{equation}

En la figura 4 se puede ver la distribuci\'{o}n triangular
con $a = 0$, $b = 1$ y $c = 3$ obtenida a partir de transformar los n\'{u}meros
generados en la secci\'{o}n \ref{generador}.

\begin{figure}[t]
\label{fig:triangular}
\centering
\includegraphics[width=3.2in]{histoTriangular}
\caption{Histograma de variables pseudoaleatorias con distribuci\'on de probabilidad triangular con par\'ametros a=0, b=1 y c=3}
\end{figure}

\section{Sistema de propulsi\'{o}n WARP}\label{simulacionpropulsor}

El sistema de propulsi\'{o}n WARP de la nave USS Enterprise consiste en dos propulsores, alimentados mediante un N\'{u}cleo WARP, donde se llevan a cabo reacciones de aniquilaci\'on materia/antimateria, moderadas por cristales de dilitio. Dichos cristales 
son procesados en una c\'amara controlada, llamada Matriz de Dilitio.

El tiempo de operaci\'on de cada propulsor es una variable aleatoria de distribuci\'on exponencial con tiempo medio de 10 d\'ias,
defini\'{e}ndose entonces las variables $X_1$ y $X_2$ de distribuci\'on exponencial con $\lambda=1/240$ horas.

El tiempo de operaci\'on entre fallos del n\'ucleo WARP es de 72 horas, pudiendo variar linealmente hasta en 12 horas, por lo que se
define una variable aleatoria $X_3$ de distribuci\'on triangular, con $a = 60$ horas, $b = 72$ horas y $c = 84$ horas.

Finalmente, las c\'amaras de dilitio funcionan de forma tal que la c\'amara principal tiene un tiempo de operaci\'on de entre 
20 y 50 horas con distribuci\'on uniforme, y la c\'amara redundante que tiene un tiempo entre fallos de 5 a 12 horas y se encuentra siempre en operaci\'on. Se definen entonces dos variables aleatorias $X_4~U[20,50]$ y $X_5~U[5, 12]$, que corresponden a los tiempos de la c\'amara principal y la redundante respectivamente.

Para generar las variables pseudoaleatorias $X_1$, $X_2$, $X_3$, $X_4$ y $X_5$ se generan las variables $u_i$ con el algoritmo de L'Ecuyer con
las semillas que se pueden ver en la tabla 1.

\begin{table}[h]
\label{seeds}
\centering
\begin{tabular}{|c|c|c|}
\hline
 $u_1$ & 154 & 226 \\
\hline
 $u_2$ & 63 & 344 \\
\hline
 $u_3$ & 8952 & 906 \\
\hline
$u_4$ & 297 & 34978 \\
\hline
$u_5$ & 239 & 90 \\
\hline
\end{tabular}
\caption{Semillas para las variables generadas $u_i$ con L'Ecuyer}
\end{table}

El tiempo de vuelo del USS Enterprise se estima mediante la siguiente variable aleatoria:
\begin{equation}
T = min\{ max\{X_1, X_2\}, X_3, max\{X_4, X_5\} \}
\end{equation}

Con 89 simulaciones de montecarlo se obtiene una
media muestral del tiempo de vuelo de $\bar{T}=34.7410$ horas, una varianza muestral de $5.9819$ horas y un desv\'io
muestral de $2.4458$ horas.
En la figura 5 se muestra el tiempo medio de vuelo en funci\'on de la cantidad de simulaciones realizadas, y en la figura 6
se puede ver el desv\'io muestral, tambi\'{e}n en funci\'on de la cantidad de simulaciones. Se puede ver que a medida que
aumenta la cantidad de simulaciones, la media se estabiliza en un valor, y el desv\'io cae como $\frac{1}{\sqrt{n}}$.

\begin{figure}[t]
\label{fig:media}
\centering
\includegraphics[width=3.2in]{tiempoVueloMedio}
\caption{Tiempo de vuelo medio en funci\'on de la cantidad de simulaciones.}
\end{figure}

\begin{figure}[t]
\label{fig:desvio}
\centering
\includegraphics[width=3.2in]{desvioTiempoVuelo}
\caption{Desv\'io muestral del tiempo de vuelo en funci\'on de la cantidad de simulaciones.}
\end{figure}

Para decidir cortar el algoritmo en 89 simulaciones se parte de la base de que
el intervalo de confianza para el estimador del tiempo medio de vuelo es el siguiente:

\begin{equation}
\label{eq:student}
\langle T \rangle - t_{n-1,1-\alpha/2}\frac{S}{\sqrt{n}}
< \langle T \rangle <
\langle T \rangle + t_{n-1,1-\alpha/2}\frac{S}{\sqrt{n}}
\end{equation}

Entonces se pide que $2 t_{n-1,1-\alpha/2}\frac{S}{\sqrt{n}}$ sea menor que el $2.5\%$ de la media que se est\'{a} estimando.
Con un nivel de significaci\'{o}n $\alpha = 0.05$, y tomando $t_{n-1,1-\alpha/2} = 1.65$, se obtiene
que esto se cumple en la iteraci\'{o}n 89.

\section{Conclusiones}\label{conclusiones}

Los resultados de los tests que se realizan al generador de n\'{u}meros 
pseudoaleatorios de L'Ecuyer muestran que es un buen generador. Esto se debe a
 que no solo no se observan hiperplanos en los gr\'{a}ficos, sino que 
tambi\'{e}n obtuvo resultados favorables en los tests $\chi^{2}$ y 
Kolmogorov-Smironov. En conclusion, el generador de L'Ecuyer nos permite 
generar una buena secuencia de n\'{u}meros pseudoaleatorios, con 
distribuci\'{o}n uniforme.
En el an\'{a}lisis del USS Enterprise, se ve que la principal limitaci\'on del 
sistema es la Matriz de Dilitio, ya que el tiempo de operaci\'on es el menor 
entre los tiempos de los componentes de la nave. La c\'amara redundante no es 
de suma ayuda ya que el tiempo entre fallos es mucho menor que el de operaci\'on del sistema principal. Adem\'as 
del hecho que el componente funciona siempre a la par del principal y no 
cuando este falla, logra que el sistema integro base su tiempo de vuelo \'unicamente 
en la Matriz de Dilitio.

\end{document}

\begin{thebibliography}%{1}

\bibitem{IEEEhowto:kopka}
Tabla de valores cr\'iticos de la distribuci\'on $\chi^2$. http://www.medcalc.be/manual/chi-square-table.php.

\bibitem{IEEhowto:kopka}
Kolmogorov-Smirnov Test. http://www.eridlc.com/onlinetextbook/appendix/table7.htm.

\end{thebibliography}l
